
clear all;
clc;
disp('test result with Z formula using mean and sigma until the current point time');
load dataAnalysisQ_G1000k.mat

INSAMPLE_YR = 2008;
INSAMPLE_QR = 4;

inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = length(SVD)-inSampleLength;

xx =zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
avgRollingTM =zeros(outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
bin=zeros(outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
PHat=zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
filename='optimization_07q4_jun9.xls';


%To calculate the rolling avg TM and rolling credit score/bins-------------

for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:,:) = durationCount(i,j,:,:);
        xx = xx + x;
        t = (i-1)*4 + j;
        
        for ii = 1 : RATE_LIST_LENGTH-1
            for jj = 1: RATE_LIST_LENGTH
                avgRollingTM (t, ii, :) = xx(ii,:)./ sum(xx(ii,:),2);
                cdf = sum( avgRollingTM(t, ii, jj:RATE_LIST_LENGTH) );
                bin(t,ii,jj)=norminv(cdf,0,1);
            end
        end
        
    end
end

disp('finish Rolling Avg Transition Matrix and Bin');


%To calculate credit cycle index Z, real and forecast----------------------

% Website tutoring upload data from xls to matlab:
% http://blinkdagger.com/matlab/matlab-using-xlsread-to-import-excel-data/

% [fitSVD] = xlsread('fittedReg_Sample98Q1_08Q4.xlsx','B6:B49');
[fcstInvPD] = xlsread('fittedReg_SampleEndTo_08Q4_09Q3_PD.xlsx','B50:B53');


Zhat = zeros (outSampleLength, 1);

inversePD = norminv(defaultFreq);

for i = 1 : outSampleLength
    
% SVDhat = SVD (1:inSampleLength + i -1, :);
% SVDhat(inSampleLength + i, :) = fcstSVD (i, :);
% ZhatVector = - zscore (SVDhat);
% Zhat(i, :) = ZhatVector(size(ZhatVector,1),:);
% save in 'optimization_07q4_jun2_3.xls';

Zhat(i,:) = - ( fcstInvPD(i,:)- mean(inversePD(1:inSampleLength+i-1,:)) )/std(inversePD(1:inSampleLength+i-1,:));
% save in 'optimization_07q4_jun2_5.xls';


end

% To minimize the distance to get estimate of w----------------------------

W = zeros (outSampleLength, 3); % first column is for metric BFS, the second colum is for metric D2
PHat_BFS = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
PHat_D2 = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
PHat_SVD = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
R = RATE_LIST_LENGTH-1;
C = RATE_LIST_LENGTH;
yrEst= INSAMPLE_YR - 1998 + 1 + floor(INSAMPLE_QR / 4);

for i = 1: outSampleLength

qrEst= mod(INSAMPLE_QR,4)+ i;
N(:,:) = durationCount(yrEst,qrEst,:,:);
P(:,:) = stdProbQ(inSampleLength+i,:,:);
X(:,:) = bin(inSampleLength,:,:);
wantedZ = Zhat(i); 

% BFS distance
obj1=@(w)difBFS(w,N,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj1,0,0.5);
W(i,1) = foundW;
PHat_BFS(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
% sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% xlswrite(filename, PHat, sheetlist);
% figure;
% subplot(outSampleLength,1,i); fplot(obj1,[-0.5 0.5],'*')
% title('Mininization Object with BFS distance',... 
%   'FontWeight','bold')


% D2 distance
obj2=@(w)difD2(w,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj2,0,0.5);
W(i,2) = foundW;
PHat_D2(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);


% SVD
obj3=@(w)difSVD(w,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj3,0,0.5);
W(i,3) = foundW;
PHat_SVD(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
xx(:,:) = PHat_SVD(i,:,:);
sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
xlswrite(filename, xx, sheetlist);
clear xx;
end




% MAE ---------------------------------------------------------------------
mae_perf = zeros(5,outSampleLength);
mae_svd_perf = zeros(5,outSampleLength);

for i = 1: outSampleLength
    P(:,:) = stdProbQ(inSampleLength+i,:,:);
    previous (:,:) = stdProbQ(inSampleLength+i-1,:,:);
    average (:,:) = avgRollingTM(inSampleLength+i-1,:,:);
    model_BFS (:,:) = PHat_BFS(i,:,:);
    model_D2 (:,:) = PHat_D2(i,:,:);
    model_SVD (:,:) = PHat_SVD(i,:,:);
    
    
    mae_perf(1,i) = mean2 (abs (P - previous));
    mae_perf(2,i) = mae (P - average);
    mae_perf(3,i) = mae (P - model_BFS);
    mae_perf(4,i) = mae (P - model_D2);
    mae_perf(5,i) = mae (P - model_SVD);    
    
    mae_svd_perf(1,i) =  abs (mean(svd(P))- mean(svd(previous)));
    mae_svd_perf(2,i) =  abs (mean(svd(P))- mean(svd(average)));
    mae_svd_perf(3,i) =  abs (mean(svd(P))- mean(svd(model_BFS)));
    mae_svd_perf(4,i) =  abs (mean(svd(P))- mean(svd(model_D2)));
    mae_svd_perf(5,i) =  abs (mean(svd(P))- mean(svd(model_SVD)));
  
 
end

% naive1=stdProbAvgQ;
% naive2=stdProbQ(inSampleLength+i,:,:);
% % n(:,:)=N(yrp,freqjp,2:10,13);
% % nn=repmat(n,1,10);
% % nt(1,1,:,:)=nn;
% % Preal=P(:,:,2:10,2:11);
% e_mod=(nt.*(Preal-PZ).^2)./(PZ.*(I-PZ));
% e_avg=(nt.*(Preal-naive1).^2)./(naive1.*(I-naive1));
% e_pre=(nt.*(Preal-naive2).^2)./(naive2.*(I-naive2));
% perf_mod = mae(e_mod);
% perf_naive1 = mae(e_avg);
% perf_naive2 = mae(e_pre);
% % 
% 
% Parameters={'foundW','Zhat', 'perf_mod', 'perf_naive1','perf_naive2'; foundW Zhat perf_mod perf_naive1 perf_naive2}';
% xlswrite(filename, Parameters,'Parameters');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Naive1 MAE of L1,L2 NSD,D1,D2   %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L1_n1=Preal-naive1;
% L2_n1=(Preal-naive1).^2;
% NSD_n1=((Preal-naive1).^2)./(naive1);
% % D1_n1=
% % D2_n1=
% 
% perf_naive1e_L1=mae(L1_n1);
% perf_naive1e_L2=sqrt(mae(L2_n1).*90)/90;
% perf_naive1e_NSD=mae(NSD_n1);
% % perf_naive1e_D1=
% % perf_naive1e_D2=
% 
% disp('fini');

 
